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BAYESIAN  COMPOTATIONS  IN  SURVIVAL  HODELS  VIA  THE  GIBBS  SAMPLER 


LYNN  KUO. 

Depeirtment  of  Statistics. 
University  of  Connecticut. 
U.  S.  A. 


ADRIAN  F.  M.  SMITH. 
Department  of  Mathematics. 
Imperial  College  London. 
United  Kingdom. 


ABSTRACT.  Survival  models  used  in  biomedical  and  reliability  contexts 
typically  involve  data  censoring,  and  may  also  Involve  constraints  in 
the  form  of  ordered  parameters.  In  addition.  Inferential  interest 
often  focuses  on  non-linear  functions  of  natural  model  parameters. 
From  a  Bayesian  statistical  analysis  perspective.  these  features 
combine  to  create  difficult  computational  problems  by  seeming  to 
require  (multi-dimensional)  numerical  integrals  over  awkwardly  defined 
regions.  This  paper  illustrates  how  these  apparent  difficulties  can  be 
overcome,  in  both  parametric'  and  non-peucametric  settings,  by  the  Gibbs 
sampler  approach  to  Bayesian  computation. 


1 .  Introduction 

For  the  Bayesiein  statistical  analysis  of  other  than  simple  stylized 
models,  the  key  tool  for  calculation  is  (multi-dimensional)  numerical 
integration;  see.  for  example.  Smith  et  al  (1987)  for  a  review  of 

available  techniques.  However,  it  is  widely  recognized  that 
considerable  numerical  sophistication  is  typically  required  in  applying 
these  techniques,  and  that  this  has  thus  far  hampered  the  development 
of  routinely  available,  user-friendly,  Bayes ieui  computational  methods. 

This  is  particularly  true  in  the  case  of  survival  models  used  in 

biomedical  and  reliability  contexts.  Here,  features  such  as  data 
censoring,  ordered  parameters,  assumed  convexity  or  concavity  of  dis¬ 
tributions,  all  conspire  to  produce  complicatedly  constrained  regions 
over  which  numerical  integrations  are  required.  Not  surprisingly,  the 
literature  therefore  contains  very  few  instances  of  fully  BayesiEui 
analyses  in  survival  contexts  (i.e. ,  presenting  full  and  accurate 
posterior  summaries,  rather  than,  say  modal  point  estimates  or  second 
derivative  uncertainty  measures). 

Recently,  however,  Gelfand  et  ai  (1991)  have  shown  that  the  Gibbs 

sampler  approach  to  Bayesian  computation  (see,  for  example,  Gelfand  and 
Smith,  1990,  and  Gelfand  et  ai,  1990)  effectively  side-steps  the 


seeming  problems  of  awkvordly  defined  integration  regions  in  triincated 

data  and  constrained  parameter  problems,  and  provides  an  easily 
implemented  computational  procedure. 

Our  purpose  in  this  paper  is  to  illustrate  the  simplicity  and  scope  of 
the  Gibbs  sampler  for  the  routine  Bayesian  analysis  of  survival  data, 
in  both  parametric  and  non-parametric  settings.  In  Section  2,  we 

briefly  review  the  Gibbs  sampler  and  its  general  stnx;ture  for 

constrained  parameter  and  censored  data  problems.  In  Section  3,  we 
provide  a  range  of  illustrations  of  how  the  methodology  proceeds  for  a 
variety  of  parametric  models  used  in  various  survival  modelling 

contexts.  In  Section  4,  we  give  a  non-parametric  illustration  of  the 
methodology. 


2.  The  Gibbs  Sampler,  Constraints  and  Censoring 

2.1  THE  GIBBS  SAMPLER 

In  what  follows,  densities  will  be  denoted  generically  by  sqiiare 
brackets,  so  that  Joint,  conditional  and  marginal  forms  for 
random  variables  U,  V  appear,  respectively,  as  [V,  V],  and  [K]. 

Meu'glnallzation  by  Integration  is  denoted  by  ll/]  *  •  [K] .  Given  a 

collection  of  random  variable^;  with  Joint  density  [Uy  U^,  we 

shall  refer  to  [U  \U  ,  n^s],  s  »  1,  2....,k,  as  the  full  conditional 
s '  r 

dens  ities. 

The  Gibbs  ssunpler  is  a  simply  described  iterative  stochastic  simulation 
scheme,  whereby  samples  drawn  from  the  full  conditional  densities  are 
used  to  provide  summaries  of  aspects  of  the  Joint  density.  Given  an 

arbitrary  set  of  starting  values,  l/^,  a  random  variate  l/J  is 

drawn  from  [t/j  . . . ,  l/^] ,  then  a  variate  is  drawn  from 
<^3 . and  so  on  until  is  drawn  from 

completes  one  iteration  of  the  sampler  and  results  in  a  generated 
vector  (l/J, . . .  ,1/').  Repeating  this  process,  after  i  iterations  we 

IK  j  j 

arrive  at  a  generated  vector  (Uy . , .  ,U^).  It  can  be  shown  that,  under 

mild  regularity  conditions  (see,  for  example,  Gema,i  and  Geman,  1984), 
as  i-»w  this  random  vector  tends  in  distribution  to  a  random  vector 

having  the  Joint  distribution  I(/j, . . ,  ,y^]. 

One  possible  procedure  for  obtaining  summaries  of  aspects  of 
{l/j, . . .  ,1/^]  of  Interest  is  therefore  the  following.  Run  ft  independent 
parallel  replications  of  the  above  saunpling  procedure,  so  that,  for  i 


Judged  to  be  sufficiently  large.  the  resulting  generated  vectors 


•V  •  J 


from  [l/j, . . .  ,1/^1 . 


,  can  be  regarded  as  an  iid  sample  of 
Standauxi  moment,  quantile  or  density 


estimation  techniques  c^Ln  then  be  employed  to  estimate  summary  features 
of  interest. 

In  the  Bayesian  Inference  context.  is  the  Joint  posterior 

density  of  unknown  model  pau'ametere  1/^,  —  ,1/^.  Univariate  marginal 

inference  summaries  for  U  .  say.  are  simply  obtained  from 

s 

(l/^j . G®*^erally,  if  marginal  inference  summaries  are  required 

for  a  specif  ied  fxmction  of  (U^, . . .  ,U^),  an  iid  saunple  of  size  M  from 
the  corresponding  marginal  density  is  immediately  available  on  substi¬ 
tuting  the  generated  vectors  J  into  the 

functional  form.  Exploration  of  bivariate  (or  higher  dimensional) 
marginal  summaries  proceeds  in  an  obvious  manner.  For  further  detail, 
and  comments  on  the  pragmatics  of  choices  of  i  and  M,  see  Gelfand  and 
Smith  (1990,  1991),  Gelfand  et  al  (1990,  1991). 

f 

The  Gibbs  sampler  thus  provides  a  simulation-based  alternative  to 

direct  numerical  integration  methods,  and  one  which  depends  only  on  our 

capacity  to  generate  random  variates  (reasonably  efficiently)  from  the 

full  conditional  densities,  ((/  \U  ,  res).  We  shall  now  look  at  this 

s'  r 

latter  issue  in  the  context  of  constrained  parameter  and  censored  data 
problems.  Our  discussion  here  will  be  kept  to  the  minimum  necesssiry  to 
give  the  reader  an  appreciation  of  how  the  Gibbs  sampler  achieves 
crucial  simplification.  For  a  much  more  complete  discussion,  see 
Gelfand  et  ai.  (1991). 

2.2  MODELS  WITH  CONSTRAINED  PARAMETERS 


Suppose  a  parametric  model  for  data  Y  involves  a  k-dimensional 

parameter  vector  0,  constrained  to  lie  In  a  subset  of  R^.  For 

simplicity  of  exposition,  we  shall  assume  here  that  does  not  depend 
on  y  (as  would  be  the  case,  for  example,  if  some  components  of  0  were 
truncation  parameters:  see  Gelfand  et  ai,  1991).  Suppose  further  that 
[y|0],  [0]  denote  the  (unconstrained)  forms  of  likelihood  and  prior,  so 
that  the  (constrained)  Joint  posterior  for  0  »  (0j,...,0^)  is  given  by 


(0lyi 


[y|0i  [01 
r  [y|0]*[o] 


Avr,  li'M  *  Copies 


' 

Diet 

A'* nil  a 

C  .  ,  .  -5  ^ 

□□ 


Proceeding  by  direct  nunerical  integration,  we  see  that  there  is  an 
immediate  problem  in  calculating  the  normalizing  constant  smd  subse¬ 
quent  problems  in  performing  (k-1) -dimensional  integrals  (over  subsets 

of  s")  to  obtain  marginal  density  forms. 

However,  consider  now  the  full  conditional  forms  required  for  the  Gibbs 
sampler.  If  >i)  denotes  the  cross-section  of  corresponding 

to  the  constraints  on  0^  Imposed  by  for  specified  values  of  Oj,  J*i, 
we  have 

[ejy.  Oj,  j*i]  «  ty|e]*[e]  ,  j*i)  . 

Moreover,  the  constraints  region  for  0^  will  typically  be  an  Interval, 
or  a  union  of  Intervals. 


It  follows  that  the  typical  random  variate  generation  task  required  for 
the  Gibbs  sampler  in  this  case,  will  simply  be  that  of  generating  from 
specified  univariate  density  shapes  tmncated  to  intervals.  This  is  a 
relatively  straightforward  task:  in  any  case,  strikingly  easier  than 
high-dimensional  numerical  integration  over  complicated  constraint 
volumes.  '  , 


2.3  CENSORED  DATA  PROBLEMS 

Suppose  a  parametric  model  for  data  Y  «  (y . y  )  involves  a 

1  n 

k-dimensional  peirameter  vector  0,  with  likelihood  defined  by 

[y|0]  =  n  [y  |0]  . 

i  ^ 

However,  suppose  that  for  J  a  m  >  1  there  exist  Vj,  Vj  such  that, 
instead  of  observing  Yj  exactly,  we  simply  observe  that  Yj  €  {Vj,  Wj), 
so  that  the  likelihood  is  actually  given  by 


B-1 

n  (y,l0] 

j-i 


n 

n 

J*m 


(We  are  here  assuming  a  simple,  fully  specified  censoring  process  for 
convenience  of  exposition.  For  a  more  general  discussion,  see  Gelfand 
et  al.  1991). 


In  this  case,  a  moment’s  reflection  reveals  that  the  full  conditional 
forms  implied  by  the  above  likelihood  combined  with  a  prior  [0]  are 
not,  in  general,  easy  forms  to  sample  from.  In  particular,  the 


integral  terms  may  have  no  closed-form  analytic  expressions,  so  that 
standard  envelope  injection  or  ratio-of-uniforms  sampling  techniques 
are  not  readily  applicable. 

However,  suppose  we  consider  Y'  «  (Y^ . Y^)  as  additional  unknowns, 

so  that  the  unknown  model  p>arameters  are  (6,  Y'),  with  the  data  given 

by  (Y*,  V,  W).  where  Y*  =  (Y, . ^  J.  V  »  (K . V)  and  1/  = 

1  B-i  *  n 

(1/  ).  Consider  now  the  full  conditionals  required  for  the  Gibbs 

ffl  n 

sampler: 

lejY*,  V,  W.  ej,  J0i,  Y']  ,  1^1 . k  , 

[Y^|Y*,  V,  V,  e,  Y^,  s*r,  s  fe  m]  ,  r  =  m . n 

Careful  examination  of  the  conditioning  variables  reveals  that  the  full 
conditionals  for  6^, ...  ,0^  reduce  to 

[0jy,  0j.  >i]  ,  1=1 . k  , 

the  forms  that  would  l«ive  obtained  in  the  uncensored  case*  Typically, 
these  forms  present  no  difficulty  for  random  variate  generation. 

For  the  full  conditionals  for  Y  ,...,Y  ,  examination  of  the  forms 

B  n 

reveals  that  these  reduce  to 

V 

r  ^ 


naunely,  the  sampling  distributions  for  the  Y^  restricted  to  the  ranges 
(Vj,  Vj).  Again,  these  typically  present  no  difficulty  for  random 
variate  generation. 

The  trick  of  treating  censored  observations  as  unknowns  in  combination 
with  the  Gibbs  saunpler  leads  to  simple  Bayesian  calculation  strategies 
in  otherwise  intractable  problems  (see,  also.  Tanner  and  Wong,  1987, 
for  a  related  manifestation  of  the  idea).  In  the  next  section,  we 
illustrate  this  concretely  by  detailing  the  forms  of  the  Gibbs  sampler 
arising  in  a  range  of  parametric  models  used  in  various  kinds  of 
survival  studies. 


3. 


Illustrations  For  Paraoietric  Survival  Models 


3. 1  ORDERED  BINOMIAL  PARAMETERS 


Consider  conditionally  indep>endent  observations  Binomial  (n^, 

0j),  1  «  ,  where  it  is  known  that  9^  s  a.  ..s  9^  and  we 

seek  to  make  inferences  about  the  9^  (or  fiinctions,  thereof,  such  as 

9^^j-9j  or  /  ®j)-  Problems  of  this  kind  arise,  for  example, 

in  reliability  development  testing  (Smith,  1977;  Fard  and  Dietrich, 
1987),  vrtiere  stages  l,...,k  correspond  to  successive  improvements  in 
reliability. 

If  the  Joint  prior  density  is  taken  to  be  proportional  to 


n 


(1-9^) 


P 


i-1 


over  the  simplex,  =  {(9j,...,9^)  0  a  a  a. . .  a  s  1},  by 

conjugacy  the  Joint  posterior  iBjY)  has  the  same  form,  with  support  S*^, 
but  with  ttj,  replaced  by  ♦  y^,  -  Y^,  resF>ectively. 


Implementation  of  the  Gibbs  sampler  is  now  seen  to  be  very  simple.  The 
full  conditionals  are  given  by 

[9j|y,  Oj.  J*l]  -  Beta  (a^+  Y^,  P^  *  -  Y^  ,  1=1 . k  . 

restricted  to  the  interval  9j_j  a  a  (9^  =  0,  9^^j  =  1),  and 

random  veu'iate  generation  is  straightforward. 

3.2  CENSORED  REGRESSION  DATA 


Schmee  and  Hahn  (1979)  modelled  log-failure  times  of  motorettes  tested 
at  four  different  temperatures  by  a  straight-line  regression  of 
log-failure  versus  transformed  temperature.  Censoring  occurred 

whenever  a  motorette  had  not  failed  at  the  end  of  the  test  period.  The 
uncensored  case  likelihood  [y|9]  is  assumed  to  derive  from  Y..  ^  a  * 

2 

PX^  *  where  ~  N(0,  <r  ),  1  =  l,,..,k,  J  =  l,...,nj,  but  the 

actual  data,  Z,  are  given  by 


y  s  W 
Uj  i 


Y  >  W 


if 


where  Is  the  total  test  time  at  temperature  corresponding  to  X^. 

To  implement  the  Gibbs  sampler,  as  indicated  in  Section  2.3,  we  include 
the  unobserved  Y.,  (i.c. ,  those  where  Y..  >  V.)  as  further  unknowns  in 

the  model,  in  addition  to  the  basic  parameters  of  interest,  a,  ^  and  <r  . 

Given  conjugate  normal  prior  forms  for  a,^  and  an  Inverse-gamma  prior 
2 

for  <r  ,  it  is  easily  verified  that  the  full  conditional  forms  for  a,  p 
2 

Euid  <r  are  straightforwardly  identified  conjugate  forms  (normal,  normal 
and  inverse-gamma,  respectively)  obtained  as  if  all  the  Y were 

precisely  observed.  The  full  conditionals  for  the  unobserved  Y. .  are 

2 

simply  N(.a  +  PX^,  <r  ),  restricted  to  the  range  Y^j  >  V^.  Again,  random 
variate  generation  from  all  these  full  conditionals  is  unproblematic. 

3.3  TRUNCATED  BIVARIATE  NORMAL  DATA 

Consider  a  bivariate  normal  process  (X^,  Y^),  i-  l,...,n,  where  some  of 

the  Y^  are  not  observed.  One  context  in  which  such  data  arises  is  in 

paired  survival  time  studies  (using  observed  logarithms  of  survival 
times),  where  observation  (y^),  of  the  second  of  the  paired  patients  is 

terminated  when  the  first  of  the  pair  dies,  so  that  Y^  is  observed  only 

if  y .  s  X  . 

1  i 

More  precisely,  we  assume  iid  pairs  (X^,  Y^  such  that  for  i  =  l,...,n. 


X  1 

■ 

r  2 

-  N  • 

*^12 

■®2  . 

9 

•  ®'l2 

2 

*"2  J 

• 

We  observe  the  pairs  (X^,  zp  with  »  y^  if  Xj  s  X^;  otherwise  we 

observe  (X^,  •),  where  *  •  indicates  that  Y^  >  X^.  Suppose  that  the 

prior  for  (9^,  e^)  is  taken  to  be  bivariate  normal  with  mean  (pj,  n^) 

and  covariance  matrix  V,  and  that  the  prior  for  the  covariance  matrix, 

£,  say  is  taken  to  be  an  inverse-Wishart,  so  that  [£  ^]  «  W  {(pR) 
p) ,  with  all  the  hyperpsirameters  V,  p  and  R  known. 

Interest  focuses  on  marginal  inferences  for  9j,  and  £,  but, 

following  Section  2.3,  unobserved  values  of  Y^  are  also  treated  as 
unknowns  in  specifying  the  Gibbs  sampler.  Defining  ■  (X^,  Y^),  T  * 


(Tj . r^),  r  =  p  +. .+  r^),  e  =  (e^,  e^)  and  n  =  (^j,  ^ 

easily  verified  that 

[e|r,  Z.  Z]  »  NiVlrT^l.  +  V)~^  f  *  n“^Z(n"^Z  ♦  V)ti,  (nZ*^  +  , 

tz"^|r.  z.  e]  *  V{(z(r  -  e){r  -  e>'  +  p/?)"^  n  +  p) 

i  ^  ^ 

and 

(YjlXj.  Z^.  6.  rl  .  «K82  -  -  «,)■  4  -  ■ 

truncated  to  (X^,  «)  if  Z^  =  *,  with  degenerate  at  Z^  otherwise. 

The  required  random  variate  generation  is  routine. 

3.4  WEI BULL  PROPORTIONAL  HAZARDS  WITH  CENSORING 

Consider  a  survival  time  model  in  which  the  hazard  function  A(t;  Z), 
for  an  individual  with  covari^te  values  Z  at  time  t,  is  given  by 

Xit;  X)  *  pt^'^  exp(ZB)  . 

where  B  =  O^,  py...,p^)'  is  a  vector  of  unknown  regression  parameters 
Bind  p  >  0  is  the  unknown  Ueibull  shape  parameter. 

If  are  explicitly  observed  survival  times  and 

censored  (T  >  t)  lifetimes,  with  Zj  denoting  covariate  values  for  the 

Jth  case,  the  likelihood  is  given  by 


Clearly,  vdiatever  the  prior  specification,  the  resulting 
(p  2) -dimensional  posterior  is  awkward  to  handle  using  standsuxi 
numerical  integration  procedures. 

However,  it  is  easily  verified  that  the  second  partial  derivatives  of 
the  log- likelihood  with  respect  to  each  of  the  p  *  2  unknown  parameters 
are  all  non-positive  (see  Dellaportas  and  Smith,  1991).  If  the  prior 

density  is  chosen  to  be  log-concave,  it  follows  that  all  the  posterior 

full  conditionals  are  log-concave.  The  import  of  this  observation  is 
that  highly  efficient  methods  exist  for  random  variate  generation  frx>m 
log-concave  densities  (see,  in  particular,  Cilks  and  Wild,  1991),  so 


that  routine,  straightforviard  Bayes iaoi  calculation  for  widely  used 
cases  of  proportional  hazards  models  is  possible  (see  Dellaportas  and 
Smith,  1991,  for  wider  exploitation  of  log-concavity). 


4.  A  Nonparameiric  Illustration 

4.  1  INTRODUCTION 

Nonparametric  Bayesian  Inference  for  the  survival  function  with  right 
censored  data  has  been  studied  by  Siisarla  and  Wan  Ryzin  (id76),  and 
Ferguson  and  Phadla  (1979).  However,  we  often  encounter  the  situation 
where  some  observations  eire  censored  from  the  left  and  some 
observations  are  censored  from  the  right  (see  Turnbull,  1974,  for 
references  to  papers  addr-^sing  doubly  censored  data  sets  from  a 
frequentlst  perspective). 

In  this  section,  ve  study  a  nonparametric  Bayesiain  approach  to  such 
problems,  which  allows  us  to  Incorporate  prior  beliefs  and  frees  us 
from  making  a  restrictive  (parametric)  model  assumption  for  the 
survival  function.  Specifically,  we  assume  that  the  distribution 

function  F  of  survival  times  has  a  prior  given  by  .ergxison’s  (1973) 
Dirichlet  process,  D(a).  The  measure  a  can  be  written  as  NF^,  where  F^ 

is  the  prior  mean  of  F  and  F^d  ~Fq)  /  (N  +  1)  is  the  prior  variance  of 

F.  The  larger  N,  the  more  strongly  the  prior  specifies  that  F  concen¬ 
trates  around  F^. 

In  the  doubly  censored  data  case,  it  is  very  difficult  to  obtain  an 

explicit  expression  for  non-paLT'eunetrlc  Bayesian  estimators  even  in  the 
form  of  the  posterior  mean.  We  shall  show,  however,  that  the  Gibbs 
SEunpler  approach,  which  augments  the  data  by  using  latent  variables 
that  decompx^se  the  number  of  the  censored  observation  into  the  possible 
numbers  of  observations  falling  into  each  interval,  provides  a 

straightforwardly  computed  numerical  solution.  As  illustrated  in 
Section  2.3,  this  augmentation  facilitates  the  specification  of 
appropriate  full  conditional  densities,  particularly  here  for  the 
survival  functions  given  the  latent  variables.  The  iterated  saimpling 
scheme  then  allows  us  to  approximate  the  posterior  distribution  of  the 
survival  function. 

4.2  THE  MODEL 

We  shall  illustrate  the  approach  using  a  model  similar  to  that  studied 
by  Turnbull  (1974),  who  proposed  a  self-consistent  algorithm  for 
computing  the  generalized  maximum  likelihood  estimators.  Here,  we  add 
the  Dirichlet  process  prior  to  the  model. 

Let  T, ,  r_ . T  denote  the  true  survival  times  of  n  individuals  that 

12  n 


could  be  observed  precisely  If  no  censoring  were  present.  The  T ^  are 

independent  and  identically  distributed  with  distribution  F\  that  is. 

Fit)  =  PiT  s  t)  for  t  fc  0.  We  consider  the  case  that  not  all  are 

observed  precisely.  For  each  i,  we  assume  that  there  ar*e  "windows"  of 
observations  and  Wj  iV^  s  vp  that  are  either  fixed  constants  or 

random  variables  independent  of  the  {T^}.  We  observe 

Xj  =  max  I  minCT^,  V^),  V^] 

Moreover,  for  each  item,  we  also  know  whether  it  is  left-censored  with 

=  V^,  or  right-censored  with  =  W^,  or  a  precisely  observed  time 

with  X.  =  T.. 

1  1 

We  assume  that  items  (or  patients)  are  examined  at  discrete  times  (for 
example,  monthly)  and  that  there  is  a  natural  discrete  time  scale  0  < 
t,  <  <...,  t  ,  with  observed  deaths  classified  into  one  of  the  m 

Id  IB 

Intervals  (0,  t,],(t,,  t_]...,(t  t  ].  Let  S,  denote  the  number  of 
1  Id  m- 1  ."I  i 

precise  observations  (=)  in  the  period  denote  the  number 

of  left-censored  (s)  entries  at  age  and  denote  the  number  of 
right-censored  (>)  entries  at  It  is  assumed  that  the  left-censored 

entries  all  occur  at  the  end  of  age  period  (t^,  The  data  ceui 

then  be  summarized  by  the  following  tabulation; 


Type  of  obs.  \  age  (0,  tj]  (tj,  ... 


(  =  ) 
(s) 

(>) 


Let  Pj  *  Pit j)  *  1  -  Fitj)  denote  the  survival  function  evaluated  at 
tj,  so  that  the  likelihood  function  is  proportional  to 


X  \ 

n  •'’j-i  -  '■y  -  '“y  ■'  ’■/  ■ 

Let  0,  “  P,  ,  -  P,  for  J  ^  1 . B  and  let  e_,  ■  P.  The  prior 

J  J-i  J  B+ 1  B 

process  specifies  that  the  distribution  of  the  0’s  is  the  Dirichlet 

distribution 


it(e) 


c 


urr  1 

n  <»j' 


where 


for  J  = 


“j'Wfo'VV'j-i”  • 

+  1,  with  ~ 


rC 


The  posterior  distribution  of  0  =  (0,,  0^,...,0  ,  0_^-)  is  known  to  be 

a  mixture  of  Dirichlet  distributions  (see  Antoniak.  1974).  In  the  next 
section,  we  show  how  the  Gibbs  sampler  side-steps  the  need  for  direct 
computation  of  this  mixture. 

4.3  APPROXIMATION  VIA  THE  GIBBS  SAMPLER 

To  employ  the  Gibbs  sampler,  we  use  the  idea  of  Section  2.3  and 
introduce  latent  variables  that  decompose  the  numbers  of  censored 
entries  into  the  numbers  of  observations  belonging  to  individual 
intervals.  Let  Z  .,  Z_,,  ...,Z  denote  the  random  variables  that  count 

the  number  of  observations  in  Pj  that  might  fall  in  the  Intervals  (0, 

^2^ . Zjj. 

Further,  let  Z,  ,,,  ...,Z _ _  ,  denote  the  number  of  observations  in  A, 

J+lf  m+lj  J 

that  might  fall  in  the  intervals  (t,,  t,^,) . it _ 

°  J  j+1  B-1  B  B 


respectively,  so  that  A 


1-/+1 


Our  objective  is  to  suffiffl2U'ize,  via  samples  generated  form  the  Gibbs 
sampler,  the  posterior  distribution  of  0  given  the  data.  The  posterior 
full  conditional  for  0  given  the  Z’s  and  the  data,  is  easily  seen  to  be 
an  up-dated  Dirichlet  distribution  depending  only  on  the  Z’s.  The 
posterior  full  conditional  for  the  Z’s  given  0  and  the  data,  is  easily 
seen  to  be  a  product  of  multinomial  distributions.  Thus,  suppose  at 
the  1th  Iteration  step  of  the  Gibbs  SBunpler,  we  have  the  realization 

0^  «  ie{.  oi, . .  ■  with  0^.  »  1.  We  then  up-date  the  Z 

1  2  ml  J 

vsu'iables  from  the  multinomial  distributions  as  follows.  For  each 
J.  J  =  l,...,m,  we  sample  . Zj*j^  from  the  multinomial  dlstrlbu- 


tion  with  saunple  size  fij  and  paLraineters 

/  5^  for  i  =  1 . J.  Similarly. 

J=1  '' 


where 


z‘^\j . 


- 

from 


the  multinomial  distribution 


with  sample  size 


and 


p^Lrameters 


•  •  •  •  '■jj  “  ®j  /  j  j  *  J  *  ^ . m+1. 

Having  sampled  the  Z  random  variables,  we  then  generate  new  0 
variables  from  the  Dirichlet  distribution  as  follows.  We  compute,  for 
each  J,  J=l,...,ii+1, 


c^l*  Si* 


and  then  sample  . 0^^^, 

1  B 

with  parameters  . 

1  c+1 


Cl> 


from  the  Dirichlet  distribution 


By  running  M  parallel  Independent  replications  of  the  sampler,  after 

the  1th  Iteration,  we  have  . and  for 

s  =  The  posterior  distribution  of  0^  for  J  =  l,...,in  +  1  can 

then  be  approximated  (for  sufficiently  large  J)  by 


f(ejdata)  -  f  Betal^.  j'  Y^)  . 

s=l  k*l 

where  Betaia.b)  denotes  the  beta  density  with  p>araffleters  a  Euid  b.  A 
posterior  estimate  of  the  0^  is  then  given  by 


0.  =  H 


rl 


H 

ii 


is 


rm^l  „1 


^1.1 


is 


Other  posterior  summaries  can  be  computed  similarly  from  the  replicated 
samples,  i  and  M  having  been  selected  to  achieve  "convergence"  to 
"smooth"  estimates. 

4.4  A  NUMERICAL  EXAMPLE 

To  illustrate  the  Gibbs  sampler  technique,  we  shall  reeuialyze  the  data 
set  given  by  Kaplan  and  Meier  (1958).  The  data  consist  of  deaths 
occurring  at  .8,  3.1,  5.4  and  9.2  months  and  losses  occurring  at  1.0, 
2.7,  7.0  amd  12.1  months.  For  comp>arison  purposes,  we  consider  the 
same  prior  sp>ecif ications  used  by  Susarla  amd  Vam  Ryzin  (1976)  in  their 


with  ^ 


Bayesian  reanalysls  of  the  data.  That  is,  =  1  "  e 

.  12  and  N  -  4,8,  and  16. 

To  apply  the  Gibbs  sampler  approau^h,  we  divide  the  positive  real  line 

into  the  following  Intervals:  (0,  .8~],  (.8~,  .8],  (.8,  1],  (1,  2.7], 

(2.7,  3.1"],  (3.r,  3.1],  (3.1,  5.4"],  (5.4',  5.4],  (5.4,  7], 

(7,  9.2'],  (9.2",  9.2],  (9.2,  12.1],  and  (12.1,  «).  We  label  these 
intervals  by  (0,  t^],  (t^,  t^] . tj^],  and  let  e^,  .  and 

Sjg,  respectively,  denote  the  probabilities  assigned  to  the  intervals. 


where  C  is  the  normalizing  constant. 


Note  that  0^,  6g,  0g,  0jj  and  in  the  likelihood  combine  simply  with 

the  corresponding  0  variables  in  the  prior  distribution,  so  that  the 
parameters  6g,  0g,  0^^  and  0jg  are  each  up-dated  by  1  in  the 

pxjsterior  distribution.  Therefore,  we  need  only  Introduce  three  Z 

variables  for  the  Incomplete  data,  namely,  Z  »  (Z„,,  Z_ .  Z,„  ,], 

1  41  ol  1J» 1 


^^52’  ^62 . ^13,2^’  ^3 


(Z  Z  Z  Z  1 

‘  10,3*  11,3'  ^12,3’  ^13,3^ 


then  sample  Zj,  for  J  =  1,  2,  and  3,  from  the  appropriate  multinomial 


distribution  with  sample  size  1  and  rescaled  probabilities. 


To  estimate  the  survival  function  at  tj,  we  accumulate  the  0j  for  i  > 

J.  For  t  between  tj  and  Interpolation  formula  that  connects 

the  survival  function  at  the  two  end  points  according  to  the  prior 
shape  can  be  used.  Tables  1  and  2  exhibit  the  Gibbs  sampler  results 
for  the  survival  function  evaluated  at  tj  with  M  «  1000  and  H  «  4000, 

both  with  i  «  10.  The  exact  Bayes  solutions  given  by  Susarla  and  Van 


Ryzin  are  also  listed  for  comparison.  The  tables  show  that  the  Gibbs 


seunpler  results  for  M  ■  1000  are  already  very  accurate  in  approxima- 


ting  the  exact  Bayes  rules.  SimilEU'  results  hold  for  N  =  16.  For 
further  illustration  of  the  Gibbs  seunpler  nethodology,  see  Kuo  (1991), 
who  reainalyses  data  from  Turnbull  (1974). 


Table  1:  Gibbs  Approximation  to  the  Bayes  Estimates  for  N  =  A 


Statistics  \  age(t) 

.8“ 

.8 

1 

2.7 

3.  l“ 

3.  1 

A 

Pj  With  M  =  1000 

.970 

.886 

.879 

.819 

.805 

.702 

with  H  =  4000 

.970 

.886 

.879 

.819 

.805 

.701 

Exact  Bayes 

.970 

.886 

.879 

.819 

.805 

.701 

Statistics  \  age(t) 

5.4” 

5.4 

7 

9.2“ 

9.2 

12.  1 

Pj  with  H  =  1000 

.632 

.529 

.491 

.437 

.305 

.253 

P^  with  «  =  4000 

.632 

.529 

.491 

.438 

.307 

.256 

Exact  Bayes 

.632 

.528 

.490 

.438 

.306 

.255 

Table  2:  Gibbs  Approximation 

to  the  Bayes 

Estimates 

for  N 

=  8 

Statistics  \  age(t) 

.8- 

.8 

1 

2.7 

3.  r 

3.  1 

P^  with  M  »  1000 

,954 

.892 

.881 

.792 

.773 

.698 

P^  with  H  »  4000 

.954 

.892 

.881 

.792 

.773 

.700 

EIxBLCt  Bayes 

.954 

.892 

.881 

.793 

.773 

.699 

Statistics  \  age(t) 

5.4“ 

5.4 

7 

9.2“ 

9.2 

12.  1 

Pj  with  M  «  1000 

.600 

.527 

.474 

.405 

.316 

.249 

Pj  with  H  »  4000 

.602 

.529 

.474 

.405 

.318 

.250 

Exact  Bayes 

.602 

.528 

.474 

.405 

.318 

.250 
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